############################################################################
# Bypassing the Enemy (2017)
# This file contains the replication code for the paper 
# "Bypassing the Enemy: Distributive Politics, Credit Claiming, and Nonstate Organizations in Brazil"
############################################################################

# P-value plots
# RDD Estimates Plot (accompanied by estimator and opt. bd functions)
# RDD Table
# P-value for t-tests 
# P-values for F-tests
# RDD Table
# External Validity Summary
# Edge function (from Green et al, 2009)
# Modified local linear plot
# RDD Balance Table 
# RD Plot
# RDD by Election Cycle
# inv_hyp function for log transformations
# rdd.table.balance.poly
# conjoint plot

#Manipulation of Running Variable Set of Functions from Cattaneo et al 2017

############################################################################

#P-value plot (Figure 1)
pv.plot <- function(x, y, h.opt=0.01, names=names(y)[1], cexpoint=1){
  plot(x, y, type="p", pch=20, cex=cexpoint, cex.main=2, axes=F, 
       ylim=c(0, 1), xlab="Vote margin (%) - Distance from Cutoff", 
       ylab="", main=names, cex.lab=1.8)
  axis(2, at=c(0.05, 0.1, 0.5, 1), labels=c("0.05", "0.1", "0.5", "1"), cex.axis=1.8, las=2)
  axis(1, at=x, labels=as.character(round(x*100, digits=2)) , cex.axis=1.4)
  abline(h=0.05, col="gray20", lty=3, lwd=3)
  abline(h=0.1, col="gray20", lty=3, lwd=3)
}

#Estimator Functions

#Robust clustered standard errors
robust.se <- function(model, cluster){
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- model$rank 
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum))
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  rcse.se <- coeftest(model, rcse.cov)
  return(list(rcse.cov, rcse.se))
}

#Difference of mean (using OLS with robust standard errors)
dom <- function(rescaled, treat, outcome, clusterT=F, cluster.var=NULL){
  model <- lm(outcome~treat)
  est <- NULL
  if (clusterT==F){
    est <- c(model$coefficients[2], sqrt(diag(vcovHC(model,type="HC3"))[2]))
    return(as.numeric(est)) 
  }
  else{
    tmp <- robust.se(model, cluster.var)[[2]]
    est <- c(tmp[2,1], tmp[2,2])
    return(est)
  }
}

#Regression with fourth-degree polynomials
polynomial <- function(rescaled, treat, outcome, deg=4, clusterT=F, cluster.var=NULL){
  model <- lm(outcome ~ treat + poly(rescaled, degree=deg,raw=T) +
                poly(rescaled*treat, degree=deg,raw=T))
  est <- NULL
  if (clusterT==F){
    est <- c(model$coefficients[2], sqrt(diag(vcovHC(model,type="HC3"))[2]))
    return(as.numeric(est))
  }
  else{
    tmp <- robust.se(model, cluster.var)[[2]]
    est <- c(tmp[2,1], tmp[2,2])
    return(est)  
  }
}

### MODIFIED CODE FROM:
### Replication code for "Testing the Accuracy of Regression Discontinuity Analysis 
### Using Experimental Benchmarks" by Donald P. Green, Terence Y. Leong, Holger L. Kern,
### Alan S. Gerber, and Christopher W. Larimer

# Local linear regression with triangular (edge) kernel

loc_lin <- function(outcome, rescaled, cutoff, running, h=windows, treat, clusterT=F,
                    cluster.var=NULL)   {
  
  temp <- (running - cutoff) / h
  kern <- as.numeric(abs(temp) <= 1) * (1 - abs(temp))
  
  runningCen <- running - cutoff
  RunningCenTreat <- I(runningCen*treat)
  sel <- kern > 0
  
  wls.pool <- lm(outcome ~ treat + runningCen + RunningCenTreat, weights = kern, subset = sel)
  est <- NULL
  
  if (clusterT==F){
  est <- c(summary(wls.pool)$coeff[2,1],sqrt(diag(vcovHC(wls.pool, type = "HC3")))[2])
  return(as.numeric(est))
  }
  else{
    cluster.var1 <- cbind(running, cluster.var)
    cluster <- cluster.var1[abs(cluster.var1[,1]) <= h,]
    tmp <- robust.se(wls.pool, cluster[,2])[[2]]
    est <- c(tmp[2,1], tmp[2,2])
    return(est)
  }
}

### CODE FROM:
### Replication code for "Testing the Accuracy of Regression Discontinuity Analysis 
### Using Experimental Benchmarks" by Donald P. Green, Terence Y. Leong, Holger L. Kern,
### Alan S. Gerber, and Christopher W. Larimer

#Optimal Bandwidth (does not account for clustering)
optimal.bw <- function(Y,X,c,reg)   {
  
  # STEP 1
  above <- X >= c
  N <- length(X)
  S2x <- sum((X-mean(X))^2) / (N-1)
  h1 <- 1.84 * sqrt(S2x) * N^(-1/5)
  
  Nh1.neg <- length(X[X < c & X >= (c-h1)])
  Nh1.pos <- length(X[X >= c & X <= (h1+c)])
  
  Yh1.neg <- sum(Y[X < c & X >= (c-h1)]) / Nh1.neg
  Yh1.pos <- sum(Y[X >= c & X <= (h1+c)]) / Nh1.pos
  
  fhatx <- (Nh1.neg + Nh1.pos) / (N*h1*2)
  
  fhatx.temp1 <- sum((Y[X < c & X >= (c-h1)] - Yh1.neg)^2)
  fhatx.temp2 <- sum((Y[X >= c & X < (h1+c)] - Yh1.pos)^2)
  
  sigmahat2 <- (fhatx.temp1 + fhatx.temp2) / (Nh1.neg + Nh1.pos)
  
  
  # STEP 2
  
  medianX.pos <- median(X[X >= c])
  medianX.neg <- median(X[X < c])
  
  keep <- !(X < medianX.neg | X > medianX.pos)
  
  X.cen <- X - c
  X.cen2 <- X.cen^2
  X.cen3 <- X.cen^3
  
  lm.out <- lm(Y ~ above + X.cen + X.cen2 + X.cen3, subset = keep)
  summary(lm.out)
  
  mhat3 <- summary(lm.out)$coef[5,1] * 6
  
  
  N.pos <- sum(above)
  N.neg <- length(above) - sum(above)
  
  h2.pos <- (sigmahat2 / (fhatx * max(mhat3^2, .01)))^(1/7) * 3.56 * N.pos^(-1/7)
  h2.neg <- (sigmahat2 / (fhatx * max(mhat3^2, .01)))^(1/7) * 3.56 * N.neg^(-1/7)
  
  
  Yhat.pos <- Y[X >= c & X <= (h2.pos+c)]
  Xhat.pos <- X[X >= c & X <= (h2.pos+c)]
  
  Yhat.neg <- Y[X < c & X >= (c-h2.neg)]
  Xhat.neg <- X[X < c & X >= (c-h2.neg)]
  
  N2.pos <- length(Xhat.pos)
  N2.neg <- length(Xhat.neg)
  
  Xhat.pos.cen <- Xhat.pos - c
  Xhat.pos.cen2 <- Xhat.pos.cen^2
  
  Xhat.neg.cen <- Xhat.neg - c
  Xhat.neg.cen2 <- Xhat.neg.cen^2
  
  
  lm.out <- lm(Yhat.pos ~ Xhat.pos.cen + Xhat.pos.cen2)
  summary(lm.out)
  mhat2.pos <- summary(lm.out)$coef[3,1] * 2
  
  lm.out <- lm(Yhat.neg ~ Xhat.neg.cen + Xhat.neg.cen2)
  summary(lm.out)
  mhat2.neg <- summary(lm.out)$coef[3,1] * 2
  
  
  # STEP 3
  
  rhat.pos <- (720 * sigmahat2) / (N2.pos * h2.pos^4) # NOT identical to IK paper
  rhat.neg <- (720 * sigmahat2) / (N2.neg * h2.neg^4) # NOT identical to IK paper
  
  if(reg == F)    { rhat.pos <- 0; rhat.neg <- 0 }    # turn off regularization
  
  
  hhat.opt <- ((2 * sigmahat2) /
                 (fhatx * ((mhat2.pos - mhat2.neg)^2 + (rhat.pos + rhat.neg))))^(1/5) *
    3.4375 * N^(-1/5)
  
  
  return(hhat.opt)
  
}

############################################################################
# Plotting function
RD_plot <- function(
  
  # Data
  running=running,
  treat=treat,
  outcome=outcome,
  cutoff=0,
  clusterT=T,
  cluster.var=cluster.var,
  min_running=F,
  max_running=F,
  
  # Intervals and Opt bw
  opt_bw=F, # opt bw
  nr_windows=50,
  
  # Estimator/s
  # Main Estimator
  main_est="dom",
  # estimator
  #confidence intervals
  ci="95%", #option: ci="90%", ci="F"
  add_est=F, # no additional estimator
  #cluster=T,
  
  #Plot aesthetics
  nr_obs=T,
  nr_obs_lab,
  label_x="Absolute distance from the cutoff",
  label_y="Estimated effect",
  plot_label="RD plot",
  label_size=1.2,
  main_size=2,
  legend=F # Keep false if desired output is grid of plots
){
  
  # WARNINGS
  if (cutoff==0) print("Assuming that cutoff is at the value of 0 for the running variable. \n 
                       To change this, input new value for the cutoff argument.")
  if (cutoff!=0) print("Rescaling the running variable to place cutoff at zero.")
  library("sandwich")
  print("This function requires the sandwich package")
  
  if(clusterT==F){
  
  # DATA
  # 0. rescaling running var and organizing data 
  rescaled <- running - cutoff
  abs_running <- abs(rescaled)
  data <- cbind(rescaled, treat, outcome, abs_running, running)
  
  # 1. get windows from intervals and add optimal bandwidth
  min_running <- ifelse(min_running==F, 0, min_running-cutoff)
  max_running <- ifelse(max_running==F, max(abs_running), max_running-cutoff)
  windows <- seq(min_running, max_running,by=(max_running/nr_windows))[-1]
  if (opt_bw!=F) windows <- sort(c(windows, opt_bw))
  if (opt_bw!=F) print("Assuming optimal bandwidth measured as distance from the cutoff.")
  
  # 2. calculate estimate(s) and CI
  # 2.a main estimate
  if (main_est=="dom") main <- dom
  if (main_est=="poly") main <- polynomial
  if (main_est=="loc_lin") main <- loc_lin
  
  # 2.b additional
  if (add_est=="dom") second <- dom
  if (add_est=="poly") second <- polynomial
  if (add_est=="loc_lin") second <- loc_lin
  
  #2.c loop
  ests <- matrix(NA,length(windows),3)
  for (i in 1:length(windows)){
    # select data
    temp <- as.data.frame(data[abs_running<=windows[i],])
    # get estimates
    if (main_est=="loc_lin"){ 
      ests[i,1:2] <- with(temp, main(rescaled=rescaled, treat=treat,
                                     outcome=outcome, running=running,
                                     cutoff=0, h=windows[i]))}
    else{
      ests[i,1:2] <- with(temp, main(rescaled=rescaled, treat=treat, outcome=outcome))}
    
    if (add_est==F) next 
    if (add_est=="loc_lin"){
      ests[i,3] <- with(temp, second(rescaled=rescaled, treat=treat,
                                     outcome=outcome, running=running,
                                     cutoff=0, h=windows[i]))[1]}
    else{    
      ests[i,3] <- with(temp, second(rescaled=rescaled, treat=treat, outcome=outcome))[1]}
    }
  }
  else{
    # DATA
    # 0. rescaling running var and organizing data 
    rescaled <- running - cutoff
    abs_running <- abs(rescaled)
    data <- cbind(rescaled, treat, outcome, abs_running, running, cluster.var)
    
    # 1. get windows from intervals and add optimal bandwidth
    min_running <- ifelse(min_running==F, 0, min_running-cutoff)
    max_running <- ifelse(max_running==F, max(abs_running), max_running-cutoff)
    windows <- seq(min_running, max_running,by=(max_running/nr_windows))[-1]
    if (opt_bw!=F) windows <- sort(c(windows, opt_bw))
    if (opt_bw!=F) print("Assuming optimal bandwidth measured as distance from the cutoff.")
    
    # 2. calculate estimate(s) and CI
    # 2.a main estimate
    if (main_est=="dom") main <- dom
    if (main_est=="poly") main <- polynomial
    if (main_est=="loc_lin") main <- loc_lin
    
    # 2.b additional
    if (add_est=="dom") second <- dom
    if (add_est=="poly") second <- polynomial
    if (add_est=="loc_lin") second <- loc_lin
    
    #2.c loop
    ests <- matrix(NA,length(windows),3)
    for (i in 1:length(windows)){
      # select data
      temp <- as.data.frame(data[abs_running<=windows[i],])
      # get estimates
      if (main_est=="loc_lin"){ 
        ests[i,1:2] <- with(temp, main(rescaled=rescaled, treat=treat, clusterT=T, cluster.var=cluster.var,
                                       outcome=outcome, running=running,
                                       cutoff=0, h=windows[i]))}
      else{
        ests[i,1:2] <- with(temp, main(rescaled=rescaled, treat=treat, outcome=outcome, cluster=T, 
                                       cluster.var=cluster.var))}
      
      if (add_est==F) next 
      if (add_est=="loc_lin"){
        ests[i,3] <- with(temp, second(rescaled=rescaled, treat=treat,
                                       outcome=outcome, running=running,clusterT=T, cluster.var=cluster.var,
                                       cutoff=0, h=windows[i]))[1]}
      else{    
        ests[i,3] <- with(temp, second(rescaled=rescaled, treat=treat, outcome=outcome, clusterT=T, 
                                       cluster.var=cluster.var))[1]}
    }  
  }
  
  #2.d confidence intervals for main estimate
  if (ci=="95%") CI <- cbind(ests[,1]+1.96*ests[,2],ests[,1]-1.96*ests[,2])
  if (ci=="90%") CI <- cbind(ests[,1]+1.64*ests[,2],ests[,1]-1.64*ests[,2])
  # 2.b Define max and min for y axis
  if (add_est!=F){
    lim_y <- NA
    lim_y[1] <- ifelse(max(CI)>=max(ests[,3]), max(CI), max(ests[,3]))
    lim_y[2] <- ifelse(min(CI)<=min(ests[,3]), min(CI), min(ests[,3]))}
  if (add_est==F){ lim_y <- c(max(CI), min(CI))}
  lim_y[3] <- lim_y[2] - (lim_y[1]-lim_y[2])/7
  lim_y[4] <- lim_y[1] - (lim_y[1]-lim_y[2])/11
  
  # 2.c Define max for x axis
  if (abs(opt_bw) > abs(max_running)) max_axis <- opt_bw
  if (abs(opt_bw) < abs(max_running)) max_axis <- max_running
  
  # 3. ordering to know nr of observations by value of running variable
  data <- as.data.frame(data[order(abs_running),])
  if (nr_obs==T) {
    if (missing(nr_obs_lab))
      stop("Missing specification for number of observations to be included in the plot. \n
           If no axis for the number of observations is needed, use nr_obs=F.")
    nr_obs_lab <- cbind(nr_obs_lab, data$abs_running[nr_obs_lab])}
  
  
  # PLOT
  # Main line (add types, mins and maxs)
  if(legend==F){
    par(mar=c(4.1,6.1,6.1,2.1))
    plot(windows, ests[,1], pch=16, axes=T, 
         #limits
         ylim=c(lim_y[3],lim_y[1]),
         xlim=c(0, max_axis),
         #labels
         xlab=label_x, ylab=label_y, main=plot_label,
         bty='n', cex.main=main_size,
         cex.lab=label_size, 
         cex=1.2)
    #Zero effect line
    abline(h=0, lty=4, lw=3)
    #opt bw vertical line
    if(opt_bw!=F){
      abline(v=opt_bw, lty=5, col="gray80", lw=3)}
    
    # Add estimates line
    lines(windows, ests[,3], lwd=3,col="gray30")
    # CIs 
    if(ci!=F){
      lines(windows, CI[,1], lty=2,col="gray30")
      lines(windows, CI[,2], lty=2,col="gray30")}
    # Number of observations
    if (nr_obs==T) {
      axis(3, at=c(nr_obs_lab[,2]), labels=c(nr_obs_lab[,1]), cex=.6, col="black", 
           lwd = 0.5, padj=1, line=1, cex.axis=.7, col.axis="black")
      mtext("Number of observations", side=3, col="black", cex=.7, adj=0)}
  }
  
  if(legend==T){
    layout(rbind(1,2), heights=c(7,1))
    par(mar=c(4.1,6.1,6.1,2.1))
    plot(windows, ests[,1], pch=16, axes=T, 
         #limits
         ylim=c(lim_y[3],lim_y[1]),
         xlim=c(0, max_axis),
         #labels
         xlab=label_x, ylab=label_y, bty='n', 
         cex.lab=label_size, 
         cex=1.2)
    # title
    title(plot_label, cex.main=main_size, line=4)
    #Zero effect line
    abline(h=0, lty=4, lw=3)
    #opt bw vertical line
    if(opt_bw!=F){
      abline(v=opt_bw, lty=5, col="gray80", lw=3)}
    
    # Add estimates line
    lines(windows, ests[,3], lwd=3,col="gray30")
    # CIs 
    if(ci!=F){
      lines(windows, CI[,1], lty=2,col="gray30")
      lines(windows, CI[,2], lty=2,col="gray30")}
    # Number of observations
    if (nr_obs==T) {
      axis(3, at=c(nr_obs_lab[,2]), labels=c(nr_obs_lab[,1]), cex=.6, col="black", 
           lwd = 0.5, padj=1, line=1, cex.axis=.7, col.axis="black")
      mtext("Number of observations", side=3, col="black", cex=.7, adj=0)}
    
    # Legend
    par(mar=c(0.5,0,0,0))
    plot.new()
    if (main_est=="dom") main_leg <- "Difference of Means"
    if (main_est=="poly") main_leg <- "Polynomial Regression"
    if (main_est=="loc_lin") main_leg <- "Local Linear Regression"
    
    if (add_est!=F){
      if (add_est=="dom") second_leg <- "Difference of Means"
      if (add_est=="poly") second_leg <- "Polynomial Regression"
      if (add_est=="loc_lin") second_leg <- "Local Linear Regression"}
    
    if (ci=="95%") ci_leg <- paste("95% Confidence Intervals ", "(", main_leg, ")", sep="")
    if (ci=="90%") ci_leg <- paste("90% Confidence Intervals ", "(", main_leg, ")", sep="")
    
    if (add_est!=F){
      legend("center", "groups", c(main_leg, second_leg, ci_leg), pch=c(16, NA, NA), 
             lty=c(NA, 1, 2), lwd=c(2, 2, 2),
             col=c("black", "gray", "gray"), bty="n", cex=0.8)}
    
    if (add_est==F){
      legend("center", "groups", c(main_leg, ci_leg), pch=c(16, NA), lty=c(NA, 2), lwd=c(2, 2),
             col=c("black", "gray"), bty="n", cex=0.8)}
  }
  
}


############################################################################
# RDD Table

rdd.table <- function(runvar = NULL, 
                      treat = NULL,
                      outcome.1 = NULL, 
                      non.zero = NULL, 
                      p = NULL, q = NULL, kernel = NULL,
                      opt.pc = NULL, clusterT = F, cluster.var = NULL){
  
  if (clusterT==F){
    #Getting difference of means
    temp.all <- as.data.frame(cbind(runvar, treat, outcome.1, non.zero))
    bw <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)    
    est <- matrix(NA,10,9)
    for (i in 1:length(bw)){
      
      temp <- as.data.frame(cbind(runvar, treat, outcome.1,
                                  non.zero))[abs(runvar)<=bw[i],]
      
      reg.1 <- with(temp, lm(outcome.1 ~ treat))
      
      est[1:3, i] <- c(summary(reg.1)$coefficients[2,1], 
                       sqrt(diag(vcovHC(reg.1, type = "HC1")))[2],
                       abs(summary(reg.1)$coefficients[2,1]/sqrt(diag(vcovHC(reg.1, type = "HC1"))))[2])
      est[4, i] <- pnorm(abs(summary(reg.1)$coefficients[2,1]/sqrt(diag(vcovHC(reg.1, type = "HC1"))))[2],
                         lower.tail = FALSE)
      est[5,i] <- with(temp, mean(outcome.1, na.rm=T))
      est[6,i] <- with(temp, sd(outcome.1, na.rm=T))
      est[7,i] <- with(temp, sd(temp$outcome.1[treat==0], na.rm=T))
      
      est[8,i]<- nrow(temp %>% filter(non.zero == 1))
      est[9,i] <- nrow(temp)
      
    } 
    
    #Kernel, optimal bd per capita
    reg.kernel.pc <- with(temp.all, edge(Y = outcome.1, X = runvar, 
                                         c = 0, h = opt.pc)) 
    est[1:3, 7] <- c(reg.kernel.pc$theta, 
                     reg.kernel.pc$robust.se[2],
                     abs(reg.kernel.pc$theta/reg.kernel.pc$robust.se[2]))
    
    temp.opt.pc <- with(temp.all, (temp.all[abs(temp.all$runvar) <= opt.pc,]))
    est[4,7]  <-  pnorm(abs(reg.kernel.pc$theta/reg.kernel.pc$robust.se[2]),
                  lower.tail = FALSE)
    est[5,7] <- with(temp.opt.pc, mean(outcome.1, na.rm=T))
    est[6,7] <- with(temp.opt.pc, sd(outcome.1, na.rm=T))
    est[7,7] <- with(temp.opt.pc, sd(temp.opt.pc$outcome.1[treat==0], na.rm=T))
      
    est[8,7] <- nrow(temp.opt.pc %>% filter(non.zero == 1))
    est[9,7] <- nrow(temp.opt.pc)
    
    #RD robust
    robust <- with(temp.all, rdrobust(y=outcome.1, x=runvar, p=p, q=q, kernel=kernel))
    
    #Bias-corrected
    
    est[1:3, 8] <- c(summary(robust)$coefficients[2,1], 
                     summary(robust)$coefficients[2,2],
                     robust$z[2])
    temp.rd <- as.data.frame(cbind(runvar, treat, outcome.1,
                non.zero))[abs(runvar)<=as.numeric(robust$tabl2.str[5,1]),]
    est[4,8] <- pnorm(abs(robust$z[2]), lower.tail=FALSE)
    est[5,8] <- with(temp.rd, mean(outcome.1, na.rm=T))
    est[6,8] <- with(temp.rd, sd(outcome.1, na.rm=T))
    est[7,8] <- with(temp.rd, sd(temp.rd$outcome.1[treat==0], na.rm=T))
    
    est[8,8] <- nrow(temp.rd %>% filter(non.zero == 1))
    est[9,8] <- nrow(temp.rd)
    
    
    #Robust  
    est[1:3, 9] <- c(summary(robust)$coefficients[3,1], 
                     summary(robust)$coefficients[3,2],
                     robust$z[3])
    temp.rd <- as.data.frame(cbind(runvar, treat, outcome.1,
               non.zero))[abs(runvar)<=as.numeric(robust$tabl2.str[5,1]),]
    est[4,9] <- pnorm(abs(robust$z[3]), lower.tail=FALSE)
    est[5,9] <- with(temp.rd, mean(outcome.1, na.rm=T))
    est[6,9] <- with(temp.rd, sd(outcome.1, na.rm=T))
    est[7,9] <- with(temp.rd, sd(temp.rd$outcome.1[treat==0], na.rm=T))
    
    est[8,9] <- nrow(temp.rd %>% filter(non.zero == 1))
    est[9,9] <- nrow(temp.rd)

  }
  else{
    #Getting difference of means
    temp.all <- as.data.frame(cbind(runvar, treat, outcome.1, non.zero, cluster.var))
    bw <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)    
    est <- matrix(NA,10,9)
    for (i in 1:length(bw)){
      
      temp <- as.data.frame(cbind(runvar, treat, outcome.1,
                                  non.zero, cluster.var))[abs(runvar)<=bw[i],]
      
      reg.1 <- with(temp, lm(outcome.1 ~ treat))
      
      est[1:3, i] <- as.numeric(robust.se(reg.1, temp$cluster.var)[[2]][2,1:3])
      est[4,i] <- pnorm(abs(as.numeric(robust.se(reg.1, temp$cluster.var)[[2]][2,3])), 
                        lower.tail=FALSE)
      est[5,i] <- with(temp, mean(outcome.1, na.rm=T))
      est[6,i] <- with(temp, sd(outcome.1, na.rm=T))
      est[7,i] <- with(temp, sd(temp$outcome.1[treat==0], na.rm=T))
      
      est[8,i]<- nrow(temp %>% filter(non.zero == 1))
      est[9,i] <- nrow(temp)
      est[10,i] <- length(unique(temp$cluster.var))
      
    } 
    
    
    #Kernel, optimal bd per capita
    reg.kernel.pc <- with(temp.all, loc_lin(outcome = outcome.1, running = runvar, 
                                            treat=treat,
                         c = 0, h = opt.pc, clusterT=T, cluster.var=cluster.var))
    est[1:3, 7] <- c(reg.kernel.pc, (reg.kernel.pc[1]/reg.kernel.pc[2]))
    est[4, 7] <- pnorm(abs(reg.kernel.pc[1]/reg.kernel.pc[2]), 
                       lower.tail=FALSE)
    temp.opt.pc <- with(temp.all, (temp.all[abs(temp.all$runvar) <= opt.pc,]))
    est[5,7] <- with(temp.opt.pc, mean(outcome.1, na.rm=T))
    est[6,7] <- with(temp.opt.pc, sd(outcome.1, na.rm=T))
    est[7,7] <- with(temp.opt.pc, sd(temp.opt.pc$outcome.1[treat==0], na.rm=T))  
    est[8,7] <- nrow(temp.opt.pc %>% filter(non.zero == 1))
    est[9, 7] <- nrow(temp.opt.pc)
    est[10, 7] <- length(unique(temp.opt.pc$cluster.var))
    
    #RD robust
    robust <- with(temp.all, rdrobust(y=outcome.1, x=runvar, p=p, q=q,
                                      kernel=kernel, cluster=cluster.var))
    
    #Bias-corrected
    
    est[1:3, 8] <- c(summary(robust)$coefficients[2,1], 
                     summary(robust)$coefficients[2,2],
                     robust$z[2])
    temp.rd <- as.data.frame(cbind(runvar, treat, outcome.1,
              non.zero, cluster.var))[abs(runvar)<=as.numeric(robust$tabl2.str[5,1]),]
    
    est[4, 8] <- pnorm(abs(robust$z[2]), lower.tail=FALSE) 
    est[5,8] <- with(temp.rd, mean(outcome.1, na.rm=T))
    est[6,8] <- with(temp.rd, sd(outcome.1, na.rm=T))
    est[7,8] <- with(temp.rd, sd(temp.rd$outcome.1[treat==0], na.rm=T))
    
    est[8,8] <- nrow(temp.rd %>% filter(non.zero == 1))
    est[9,8] <- nrow(temp.rd)
    est[10, 8] <- length(unique(temp.rd$cluster.var))

    #Robust  
    est[1:3, 9] <- c(summary(robust)$coefficients[3,1], 
                     summary(robust)$coefficients[3,2],
                     robust$z[3])
    temp.rd <- as.data.frame(cbind(runvar, treat, outcome.1,
               non.zero, cluster.var))[abs(runvar)<=as.numeric(robust$tabl2.str[5,1]),]
    
    est[4,9] <- pnorm(abs(robust$z[3]), lower.tail=FALSE)
    est[5,9] <- with(temp.rd, mean(outcome.1, na.rm=T))
    est[6,9] <- with(temp.rd, sd(outcome.1, na.rm=T))
    est[7,9] <- with(temp.rd, sd(temp.rd$outcome.1[treat==0], na.rm=T))
    
    est[8,9] <- nrow(temp.rd %>% filter(non.zero == 1))
    est[9,9] <- nrow(temp.rd)
    est[10, 9] <- length(unique(temp.rd$cluster.var))
     
  }
  
   
  colnames(est) <- c(bw*100, "Local Linear", "Local Poly Bias", "Local Poly Robust")
  rownames(est) <- c("Estimates", 
                     "Std. error", 
                     "t-stat",
                     "one-tailed p-value",
                     "Mean dep. var.", 
                     "Standard-deviation dep. var.",
                     "Standard-deviation dep. var. (control)",
                     "non zero", "n", "clusters") 
  table <- xtable(est,
                  align="r|ccccccccc", 
                  include.rownames=TRUE)
  return(table) 
}

############################################################################
# P-value  for t-tests

#function to get p-value of t-test
t.test.p <- function(y, x){
  out <- t.test(y~x) 
  return(as.numeric(out[3]))
}

############################################################################
# P-values for F-test

## function to get f-stat and p-value of f-test
f.test<- function (y=treat, x) {
  temp <- lm(y~x)
  f <- as.vector(summary(temp)$fstatistic)
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  # attributes(p) <- NULL
  final <- p
  return(final)
}

############################################################################
# EXTERNAL VALIDITY SUMMARY 

summary.stats <- function(data=data){
  n_obs <- sapply(data, function(x) sum(!is.na(x))) 
  min <- as.matrix(apply(data, 2, FUN=min, na.rm=T))
  mean <- as.matrix(apply(data, 2, FUN=mean, na.rm=T))
  median <- as.matrix(apply(data, 2, FUN=median, na.rm=T))
  max <- as.matrix(apply(data, 2, FUN=max, na.rm=T))
  sd <- as.matrix(apply(data, 2, FUN=sd, na.rm=T))
  missing <- sapply(data, function(x) sum(is.na(x)))  
  summary <- as.data.frame(cbind(n_obs, min, mean, median, sd, max, missing))
  names(summary) <- c("n obs", "min", "mean", "median", "sd", "max", "missing")
  return(summary)
}


############################################################################
### Replication code for "Testing the Accuracy of Regression Discontinuity Analysis 
### Using Experimental Benchmarks" by Donald P. Green, Terence Y. Leong, Holger L. Kern,
### Alan S. Gerber, and Christopher W. Larimer

edge <- function(Y,X,c,h)   {
  
  require(sandwich)
  
  temp <- (X - c) / h
  kern <- as.numeric(abs(temp) <= 1) * (1 - abs(temp))
  
  above <- X >= c
  X.cen <- X - c
  
  X.cenAbove <- X.cen * above
  sel <- kern > 0
  
  
  wls.pool <- lm(Y ~ above + X.cen + X.cenAbove, weights = kern, subset = sel)
  summary(wls.pool)
  
  theta <- summary(wls.pool)$coeff[2,1]
  se <- sqrt(diag(vcovHC(wls.pool, type = "HC1")))[2]
  N.unweighted <- length(wls.pool$weights)
  N.weighted <- sum(wls.pool$weights)
  N.above <- table(above[sel])[2]
  N.below <- table(above[sel])[1]
  se_fit <- predict(wls.pool, se.fit=T)$se.fit
  
  
  return(list(theta = theta, se = se, N.unweighted = N.unweighted,
              N.weighted = N.weighted, N.above = N.above, N.below = N.below,
              coeff = summary(wls.pool)$coeff, fitted = wls.pool$fitted.values, se_fit = se_fit,
              robust.se = sqrt(diag(vcovHC(wls.pool, type = "HC1")))))
  
}


############################################################################
### MODIFIED CODE FROM:
### Replication code for "Testing the Accuracy of Regression Discontinuity Analysis 
### Using Experimental Benchmarks" by Donald P. Green, Terence Y. Leong, Holger L. Kern,
### Alan S. Gerber, and Christopher W. Larimer

# Function to estimate local average treatment effects using local linear regression with
# triangular (edge) kernel

local.linear.1 <-function(outcome = outcome,  runvar = runvar , h = H,
                        ylim=NULL, ylab, xlab="Vote Margin (%)", axes=TRUE, main=main,
                        opt.bw1=0.1456851, opt.bw2=0.1731197){
   
      opt.bwcov <- optimal.bw(Y = outcome, 
                              X = runvar, 
                              c = 0, reg = T)
      est <- matrix(NA,length(h),2)
        
      for (i in 1:length(h)) {
        temp <- edge(Y = outcome, 
                     X = runvar, c = 0, h = h[i]) 
        
        est[i,1] <- temp$theta
        est[i,2] <- temp$se
        
      } 
      
      
      upper.95 <- est[,1] + (qnorm(.975)*est[,2])
      lower.95 <- est[,1] - (qnorm(.975)*est[,2])
      upper.90 <- est[,1] + (qnorm(.950)*est[,2])
      lower.90 <- est[,1] - (qnorm(.950)*est[,2])
      
      if(is.null(ylim)) {
        ylim <- c(ifelse(min(lower.95)>0,0, min(lower.95)),
                  ifelse(max(upper.95)<0,0, max(upper.95)))
      }
      
      plot(h, est[,1], type="o", pch=1, col="black", main=main,
           xlab=xlab, ylab=ylab, ylim=ylim, axes=axes, col.main="grey40",
           col.lab="gray40", col.sub="gray40", fg="gray40", lwd=2)
      
      
      if(axes==FALSE){
        axis(1, at=h, labels=as.character(H*100))
        axis(2, at=seq(min(lower.95), max(upper.95), by=(max(upper.95) - min(lower.95)/10)), 
             labels=
               as.character(seq(min(lower.95), max(upper.95), 
                                by=(max(upper.95) - min(lower.95)/10))), las=2)
      }
      
      #95%
      lines(h, 
            upper.95, # diff in means upper bound
            lty = 3, lwd = 1.5)
      #90%
      lines(h, 
            upper.90, # diff in means upper bound
            lty = 4, lwd = 1.5)
      #95%
      lines(h,
            lower.95, # diff in means lower bound
            lty = 3, lwd = 1.5)
      #90%
      lines(h, 
            lower.90, # diff in means upper bound
            lty = 4, lwd = 1.5)
      
      abline(h=0, lty=2, lwd=2, col="red")
      abline(v = opt.bw1, col = "gray", lty=5, lwd = 1.3) 
      abline(v = opt.bw2, col = "gray", lty=6, lwd = 1.3) 
      abline(v = opt.bwcov, col = "gray", lty=7, lwd = 1.3)
}
  
local.linear <-function(outcome = outcome,  runvar = runvar , h = h,
                          ylim=NULL, ylab, xlab="Vote Margin (%)", axes=TRUE){
  
  est <- matrix(NA,length(h),2)
  
  for (i in 1:length(h)) {
    temp <- edge(Y = outcome, 
                 X = runvar, c = 0, h = h[i]) 
    
    est[i,1] <- temp$theta
    est[i,2] <- temp$se
    
  } 
  
  
  upper.95 <- est[,1] + (qnorm(.975)*est[,2])
  lower.95 <- est[,1] - (qnorm(.975)*est[,2])
  upper.90 <- est[,1] + (qnorm(.950)*est[,2])
  lower.90 <- est[,1] - (qnorm(.950)*est[,2])
  
  if(is.null(ylim)) {
    ylim <- c(ifelse(min(lower.95)>0,0, min(lower.95)),
              ifelse(max(upper.95)<0,0, max(upper.95)))
  }
  
  opt <- optimal.bw(Y = outcome, X = runvar, 
                                c = 0, reg = T)
  
  plot(h, est[,1], type="o", pch=1, col="black", 
       xlab=xlab, ylab=ylab, main="", ylim=ylim, axes=axes, col.main="grey40",
       col.lab="gray40", col.sub="gray40", fg="gray40", lwd=2)
  
  
  if(axes==FALSE){
    axis(1, at=h, labels=as.character(h*100))
    axis(2, at=seq(min(lower.95), max(upper.95), by=(max(upper.95) - min(lower.95)/10)), 
         labels=
           as.character(seq(min(lower.95), max(upper.95), 
                            by=(max(upper.95) - min(lower.95)/10))), las=2)
  }
  
  #95%
  lines(h, 
        upper.95, # diff in means upper bound
        lty = 3, lwd = 1.5)
  #90%
  lines(h, 
        upper.90, # diff in means upper bound
        lty = 4, lwd = 1.5)
  #95%
  lines(h,
        lower.95, # diff in means lower bound
        lty = 3, lwd = 1.5)
  #90%
  lines(h, 
        lower.90, # diff in means upper bound
        lty = 4, lwd = 1.5)
  
  abline(h=0, lty=2, lwd=2, col="red")
  abline(v = opt, col = "gray", lwd = 1.3) 
}



Fbalance <- function(bw=bw, data=data_200315nspv2){
  data.bw <- data[abs(data$vote_margin_share) <= bw,]
  

  # 0. Select prognostic covariates 
  covariates.temp <- with(data.bw, 
                          cbind(X2000_incomepercapita.2, X2000_doctor.2, 
                                X2000_idh_education.2, X2000_idh_income.2,
                                X2000_idh_longevity.2,
                                X2000_illiterate.2, X2000_infant.2, 
                                X2000_pop.2, X2000_poverty.2, 
                                p_13.2, p_45.2, f.nom_13.2, f.nom_45.2, g_13.2, g_45.2, 
                                e.nom_13.2, e.nom_45.2,
                                Comparecimento1t.2, sex.2, age.2))
  
  p.val <- f.test(data.bw$treat, covariates.temp)  # observed p-val of f test
  
  return(p.val)}


############## Function adapted from rdplot in rdrobust package
############## Added commands for publication plots
############## (adapted by Guadalupe Tuñón)

rd.plot <- function (y, x, subset = NULL, c = 0, p = 4, numbinl = NULL, 
                     numbinr = NULL, binselect = "esmv", lowerend = NULL, upperend = NULL, 
                     scale = 1, scalel = 1, scaler = 1, hide = FALSE, cisze=.10,
                     title = NULL, x.label = NULL, y.label = NULL, x.lim = NULL, 
                     y.lim = NULL, col.dots = NULL, col.lines = NULL, type.dots = NULL, wd.lines = 3, 
                     N=TRUE, ...) 
{
  call <- match.call()
  if (!is.null(subset)) {
    x <- x[subset]
    y <- y[subset]
  }
  na.ok <- complete.cases(x) & complete.cases(y)
  x <- x[na.ok]
  y <- y[na.ok]
  if (is.null(lowerend)) {
    lowerend = min(x)
  }
  if (is.null(upperend)) {
    upperend = max(x)
  }
  x_low = lowerend
  x_upp = upperend
  if (is.null(col.lines)) {
    col.lines = "black"
  }
  if (is.null(col.dots)) {
    col.dots = 1
  }
  if (is.null(type.dots)) {
    type.dots = 20
  }
  size = sum(x >= x_low & x <= x_upp)
  y = y[x >= x_low & x <= x_upp]
  x = x[x >= x_low & x <= x_upp]
  x_l = x[x < c]
  x_r = x[x >= c]
  y_l = y[x < c]
  y_r = y[x >= c]
  x.min = min(x)
  x.max = max(x)
  range_l = max(x_l) - min(x_l)
  n_l = length(x_l)
  range_r = max(x_r) - min(x_r)
  n_r = length(x_r)
  n = n_l + n_r
  meth = "es"
  exit = 0
  if (c <= x.min | c >= x.max) {
    print("c should be set within the range of x")
    exit = 1
  }
  if (p <= 0) {
    print("p should be a positive number")
    exit = 1
  }
  if (scale <= 0 | scalel <= 0 | scaler <= 0) {
    print("scale should be a positive number")
    exit = 1
  }
  p_ceiling = ceiling(p)/p
  if (p_ceiling != 1) {
    print("p should be an integer number")
    exit = 1
  }
  if (exit > 0) {
    stop()
  }
  p1 = p + 1
  compute = 0
  rp_l = matrix(NA, n_l, p + 1)
  rp_r = matrix(NA, n_r, p + 1)
  for (j in 1:p1) {
    rp_l[, j] = x_l^(j - 1)
    rp_r[, j] = x_r^(j - 1)
  }
  gamma_p1_l = lm(y_l ~ rp_l - 1)$coeff
  gamma_p1_r = lm(y_r ~ rp_r - 1)$coeff
  mu0_p1_l = rp_l %*% gamma_p1_l
  mu0_p1_r = rp_r %*% gamma_p1_r
  J_star_orig = c(numbinl, numbinr)
  y_l.sq = y_l^2
  y_r.sq = y_r^2
  gamma_p2_l = lm(y_l.sq ~ rp_l - 1)$coeff
  gamma_p2_r = lm(y_r.sq ~ rp_r - 1)$coeff
  drp_l = matrix(NA, n_l, p)
  drp_r = matrix(NA, n_r, p)
  for (j in 1:p) {
    drp_l[, j] = j * x_l^(j - 1)
    drp_r[, j] = j * x_r^(j - 1)
  }
  mu1_hat_l = drp_l %*% (gamma_p1_l[2:p1])
  mu1_hat_r = drp_r %*% (gamma_p1_r[2:p1])
  ind.l = order(x_l)
  ind.r = order(x_r)
  x.i.l = x_l[ind.l]
  x.i.r = x_r[ind.r]
  y.i.l = y_l[ind.l]
  y.i.r = y_r[ind.r]
  dxi.l = (x.i.l[2:length(x.i.l)] - x.i.l[1:(length(x.i.l) - 
                                               1)])
  dxi.r = (x.i.r[2:length(x.i.r)] - x.i.r[1:(length(x.i.r) - 
                                               1)])
  dyi.l = (y.i.l[2:length(y.i.l)] - y.i.l[1:(length(y.i.l) - 
                                               1)])
  dyi.r = (y.i.r[2:length(y.i.r)] - y.i.r[1:(length(y.i.r) - 
                                               1)])
  x.bar.i.l = (x.i.l[2:length(x.i.l)] + x.i.l[1:(length(x.i.l) - 
                                                   1)])/2
  x.bar.i.r = (x.i.r[2:length(x.i.r)] + x.i.r[1:(length(x.i.r) - 
                                                   1)])/2
  rp.i_l = matrix(NA, n_l - 1, p + 1)
  rp.i_r = matrix(NA, n_r - 1, p + 1)
  drp.i_l = matrix(NA, n_l - 1, p)
  drp.i_r = matrix(NA, n_r - 1, p)
  for (j in 1:p1) {
    rp.i_l[, j] = x.bar.i.l^(j - 1)
    rp.i_r[, j] = x.bar.i.r^(j - 1)
  }
  for (j in 1:p) {
    drp.i_l[, j] = j * x.bar.i.l^(j - 1)
    drp.i_r[, j] = j * x.bar.i.r^(j - 1)
  }
  mu0.i_hat_l = rp.i_l %*% gamma_p1_l
  mu0.i_hat_r = rp.i_r %*% gamma_p1_r
  mu2.i_hat_l = rp.i_l %*% gamma_p2_l
  mu2.i_hat_r = rp.i_r %*% gamma_p2_r
  mu0_hat_l = rp_l %*% gamma_p1_l
  mu0_hat_r = rp_r %*% gamma_p1_r
  mu2_hat_l = rp_l %*% gamma_p2_l
  mu2_hat_r = rp_r %*% gamma_p2_r
  mu1.i_hat_l = drp.i_l %*% (gamma_p1_l[2:p1])
  mu1.i_hat_r = drp.i_r %*% (gamma_p1_r[2:p1])
  sigma2_hat_l.bar = mu2.i_hat_l - mu0.i_hat_l^2
  sigma2_hat_r.bar = mu2.i_hat_r - mu0.i_hat_r^2
  sigma2_hat_l = mu2_hat_l - mu0_hat_l^2
  sigma2_hat_r = mu2_hat_r - mu0_hat_r^2
  J.fun = function(B, V) {
    ceiling((((2 * B)/V) * n)^(1/3))
  }
  var.y_l = var(y_l)
  var.y_r = var(y_r)
  B.es.hat.dw = c(((c - x.min)^2/(12 * n)) * sum(mu1_hat_l^2), 
                  ((x.max - c)^2/(12 * n)) * sum(mu1_hat_r^2))
  V.es.hat.dw = c((0.5/(c - x.min)) * sum(dxi.l * dyi.l^2), 
                  (0.5/(x.max - c)) * sum(dxi.r * dyi.r^2))
  V.es.chk.dw = c((1/(c - x.min)) * sum(dxi.l * sigma2_hat_l.bar), 
                  (1/(x.max - c)) * sum(dxi.r * sigma2_hat_r.bar))
  J.es.hat.dw = J.fun(B.es.hat.dw, V.es.hat.dw)
  J.es.chk.dw = J.fun(B.es.hat.dw, V.es.chk.dw)
  B.qs.hat.dw = c((n_l^2/(24 * n)) * sum(dxi.l^2 * mu1.i_hat_l^2), 
                  (n_r^2/(24 * n)) * sum(dxi.r^2 * mu1.i_hat_r^2))
  V.qs.hat.dw = c((1/(2 * n_l)) * sum(dyi.l^2), (1/(2 * n_r)) * 
                    sum(dyi.r^2))
  V.qs.chk.dw = c((1/n_l) * sum(sigma2_hat_l), (1/n_r) * sum(sigma2_hat_r))
  J.qs.hat.dw = J.fun(B.qs.hat.dw, V.qs.hat.dw)
  J.qs.chk.dw = J.fun(B.qs.hat.dw, V.qs.chk.dw)
  J.es.hat.mv = c(ceiling((var.y_l/V.es.hat.dw[1]) * (n/log(n)^2)), 
                  ceiling((var.y_r/V.es.hat.dw[2]) * (n/log(n)^2)))
  J.es.chk.mv = c(ceiling((var.y_l/V.es.chk.dw[1]) * (n/log(n)^2)), 
                  ceiling((var.y_r/V.es.chk.dw[2]) * (n/log(n)^2)))
  J.qs.hat.mv = c(ceiling((var.y_l/V.qs.hat.dw[1]) * (n/log(n)^2)), 
                  ceiling((var.y_r/V.qs.hat.dw[2]) * (n/log(n)^2)))
  J.qs.chk.mv = c(ceiling((var.y_l/V.qs.chk.dw[1]) * (n/log(n)^2)), 
                  ceiling((var.y_r/V.qs.chk.dw[2]) * (n/log(n)^2)))
  if (binselect == "es") {
    J_star_orig = J.es.hat.dw
    meth = "es"
    binselect_type = "IMSE-optimal evenly-spaced method using spacings estimators"
    J_IMSE = J.es.hat.dw
    J_MV = J.es.hat.mv
  }
  if (binselect == "espr") {
    J_star_orig = J.es.chk.dw
    meth = "es"
    binselect_type = "IMSE-optimal evenly-spaced method using polynomial regression"
    J_IMSE = J.es.chk.dw
    J_MV = J.es.chk.mv
  }
  if (binselect == "esmv") {
    J_star_orig = J.es.hat.mv
    meth = "es"
    binselect_type = "mimicking variance evenly-spaced method using spacings estimators"
    J_IMSE = J.es.hat.dw
    J_MV = J.es.hat.mv
  }
  if (binselect == "esmvpr") {
    J_star_orig = J.es.chk.mv
    meth = "es"
    binselect_type = "mimicking variance evenly-spaced method using polynomial regression"
    J_IMSE = J.es.chk.dw
    J_MV = J.es.chk.mv
  }
  if (binselect == "qs") {
    J_star_orig = J.qs.hat.dw
    meth = "qs"
    binselect_type = "IMSE-optimal quantile-spaced method using spacings estimators"
    J_IMSE = J.qs.hat.dw
    J_MV = J.qs.hat.mv
  }
  if (binselect == "qspr") {
    J_star_orig = J.qs.chk.dw
    meth = "qs"
    binselect_type = "IMSE-optimal quantile-spaced method using polynomial regression"
    J_IMSE = J.qs.chk.dw
    J_MV = J.qs.chk.mv
  }
  if (binselect == "qsmv") {
    J_star_orig = J.qs.hat.mv
    meth = "qs"
    binselect_type = "mimicking variance quantile-spaced method using spacings estimators"
    J_IMSE = J.qs.hat.dw
    J_MV = J.qs.hat.mv
  }
  if (binselect == "qsmvpr") {
    J_star_orig = J.qs.chk.mv
    meth = "qs"
    binselect_type = "mimicking variance quantile-spaced method using polynomial regression"
    J_IMSE = J.qs.chk.dw
    J_MV = J.qs.chk.mv
  }
  if (scale > 1 & scalel == 1 & scaler == 1) {
    scalel = scaler = scale
  }
  J_star_l = scalel * J_star_orig[1]
  J_star_r = scaler * J_star_orig[2]
  if (!is.null(numbinl) & !is.null(numbinr)) {
    J_star_l = numbinl
    J_star_r = numbinr
    binselect_type = "manually evenly spaced"
  }
  scale_l = J_star_l/J_IMSE[1]
  scale_r = J_star_r/J_IMSE[2]
  bin_x_l = rep(0, length(x_l))
  bin_x_r = rep(0, length(x_r))
  jump_l = range_l/J_star_l
  jump_r = range_r/J_star_r
  if (meth == "es") {
    jumps_l = seq(min(x_l), max(x_l), jump_l)
    jumps_r = seq(min(x_r), max(x_r), jump_r)
  }
  if (meth == "qs") {
    jumps_l = quantile(x_l, probs = seq(0, 1, 1/J_star_l))
    jumps_r = quantile(x_r, probs = seq(0, 1, 1/J_star_r))
  }
  for (k in 1:(J_star_l - 1)) {
    bin_x_l[x_l >= jumps_l[k] & x_l < jumps_l[k + 1]] = -J_star_l + 
      k - 1
  }
  bin_x_l[x_l >= jumps_l[(J_star_l)]] = -1
  for (k in 1:(J_star_r - 1)) {
    bin_x_r[x_r >= jumps_r[k] & x_r < jumps_r[k + 1]] = k
  }
  bin_x_r[x_r >= jumps_r[(J_star_r)]] = J_star_r
  bin_xlmean = bin_ylmean = rep(0, J_star_l)
  bin_xrmean = bin_yrmean = rep(0, J_star_r)
  bin_nl = rep(0, J_star_l)
  bin_nr = rep(0, J_star_r)
  
  for (k in 1:(J_star_l)) {
    bin_xlmean[k] = mean(c(jumps_l[k], jumps_l[k + 1]))
    bin_ylmean[J_star_l - k + 1] = mean(y_l[bin_x_l == -k])
    bin_nl[J_star_l - k + 1] <- length(y_l[bin_x_l == -k])
  }
  for (k in 1:(J_star_r)) {
    bin_xrmean[k] = mean(c(jumps_r[k], jumps_r[k + 1]))
    bin_yrmean[k] = mean(y_r[bin_x_r == k])
    bin_nr[k] <- length(y_r[bin_x_r == k])
  }
  bin_x = c(bin_x_l, bin_x_r)
  bin_xmean = c(bin_xlmean, bin_xrmean)
  bin_ymean = c(bin_ylmean, bin_yrmean)
  x_sup = c(x_l, x_r)
  if (hide == "FALSE") {
    if (is.null(title)) {
      title = "RD Plot"
    }
    if (is.null(x.label)) {
      x.label = "X axis"
    }
    if (is.null(y.label)) {
      y.label = "Y axis"
    }
    if (is.null(x.lim)) {
      x.lim = c(min(x_l), max(x_r))
    }
    if (is.null(y.lim)) {
      y.lim = c(min(c(y_l, y_r)), max(c(y_l, y_r)))
    }
    
    bin_x <- bin_xmean[order(bin_xmean)]
    bin_y <- bin_ymean[order(bin_xmean)]
    bin_xr <- bin_x[bin_x >= 0]
    bin_xl <- bin_x[bin_x < 0]
    bin_yr <- bin_y[bin_x >= 0]
    bin_yl <- bin_y[bin_x < 0]
    
    wl <- sqrt(bin_nl/pi)
    wr <- sqrt(bin_nr/pi)
    wl[wl==0] <- 1
    wr[wr==0] <- 1
    
    plot(bin_x, bin_y, 
         main = title, xlab = x.label, ylab = y.label, ylim = y.lim, cex.lab=1.5, cex=1.3,
         xlim = x.lim, col = col.dots, pch = type.dots, type = "n")
    symbols(bin_xr, bin_yr, circles = wr, lwd=1.5,
            inches = cisze, fg = "navy", bg = alpha("navy", .7), add = T)
    symbols(bin_xl, bin_yl, circles = wl, lwd=1.5,
            inches = cisze, fg = "darkred", bg = alpha("darkred", .4), add = T)
    lines(x_l[order(x_l)], mu0_p1_l[order(x_l)], type = "l", 
          col = col.lines, lwd = wd.lines)
    lines(x_r[order(x_r)], mu0_p1_r[order(x_r)], type = "l", 
          col = col.lines, lwd = wd.lines)
    abline(v = c, lty=2, lwd=1.2)
    
    if (N==TRUE){
      legend(
        "bottomright", 
        legend=paste0('N=',n), 
        bty="n",
        #pch=21,
        #pt.cex=cex_brks2,
        cex = 1.2
      )
    }
    
    
  }
  tabl1.str = matrix(NA, 14, 2)
  tabl1.str[1, ] = formatC(c(n_l, n_r), digits = 0, format = "f")
  tabl1.str[2, ] = formatC(c(p, p), digits = 0, format = "f")
  tabl1.str[3, ] = formatC(c(scalel, scaler), digits = 0, format = "f")
  tabl1.str[4, ] = c("", "")
  tabl1.str[5, ] = formatC(c(J_star_l, J_star_r), digits = 0, 
                           format = "f")
  tabl1.str[6, ] = formatC(c(jump_l, jump_r), digits = 4, format = "f")
  tabl1.str[7, ] = c("", "")
  tabl1.str[8, ] = formatC(c(J_IMSE), digits = 0, format = "f")
  tabl1.str[9, ] = formatC(c(J_MV), digits = 0, format = "f")
  tabl1.str[10, ] = c("", "")
  tabl1.str[11, ] = c("", "")
  tabl1.str[12, ] = formatC(c(scale_l, scale_r), digits = 4, 
                            format = "f")
  tabl1.str[13, ] = formatC(c(1/(1 + scale_l^3), 1/(1 + scale_r^3)), 
                            digits = 4, format = "f")
  tabl1.str[14, ] = formatC(c(scale_l^3/(1 + scale_l^3), scale_r^3/(1 + 
                                                                      scale_r^3)), digits = 4, format = "f")
  rownames(tabl1.str) = c("Number of Obs.", "Polynomial Order", 
                          "Scale", "", "Selected Bins", "Bin Length", "", "IMSE-optimal bins", 
                          "Mimicking Variance bins", "", "Relative to IMSE-optimal:", 
                          "Implied scale", "WIMSE variance weight", "WIMSE bias weight")
  colnames(tabl1.str) = c("Left", "Right")
  results = matrix(NA, 10, 2)
  results[1, ] = c(n_l, n_r)
  results[2, ] = c(p, p)
  results[3, ] = c(scalel, scaler)
  results[4, ] = c(J_star_l, J_star_r)
  results[5, ] = c(jump_l, jump_r)
  results[6, ] = J_IMSE
  results[7, ] = J_MV
  results[8, ] = c(scale_l, scale_r)
  results[9, ] = c(1/(1 + scale_l^3), 1/(1 + scale_r^3))
  results[10, ] = c(scale_l^3/(1 + scale_l^3), scale_r^3/(1 + scale_r^3))
  rownames(results) = c("Number of Obs.", "Polynomial Order", 
                        "Chosen Scale", "Selected bins", "Bin Length", "IMSE-optimal bins", 
                        "Mimicking Variance bins", "Implied scale", "WIMSE variance weight", 
                        "WIMSE bias weight")
  colnames(results) = c("Left", "Right")
  coef = matrix(NA, p + 1, 2)
  coef[, 1] = c(gamma_p1_l)
  coef[, 2] = c(gamma_p1_r)
  colnames(coef) = c("Left", "Right")
  out = list(method = binselect_type, results = results, coef = coef, 
             tabl1.str = tabl1.str)
  out$call <- match.call()
  class(out) <- "rdplot"
  return(invisible(out))
}

########## RDD election cycle
rdd.table.cycle <- function(runvar=NULL, 
                      treat=NULL,
                      outcome.1= NULL, 
                      non.zero = NULL, yeare=NULL,
                      p=NULL, q=NULL, kernel=NULL){
  
  #Getting difference of means
  temp.all <- as.data.frame(cbind(runvar, treat, outcome.1, non.zero, yeare))
  est <- matrix(NA, 8, 4)
  two.sided.pvalue <- matrix(NA, 4, 1)
  one.sided.pvalue <- matrix(NA, 4, 1)
  rownames(two.sided.pvalue) <- c("Local Poly Robust, 2000", "Local Poly Robust, 2004", 
                                  "Local Poly Robust, 2008", "Local Poly Robust, 2012")
  rownames(one.sided.pvalue) <- c("Local Poly Robust, 2000", "Local Poly Robust, 2004", 
                                  "Local Poly Robust, 2008", "Local Poly Robust, 2012")
  
  #RD robust
  year <- c(2000, 2004, 2008, 2012)
  for (i in 1:length(year)){
    temp <-  temp.all[temp.all$yeare==year[i],]
    
    robust <- with(temp, rdrobust(y=outcome.1, x=runvar, p=p, q=q, kernel=kernel))
    
    #Robust  
    est[1:3, i] <- c(summary(robust)$coefficients[3,1], 
                     summary(robust)$coefficients[3,2],
                     robust$z[3])
    temp.rd <- as.data.frame(cbind(runvar, treat, outcome.1,
                                   non.zero))[abs(runvar)<=as.numeric(robust$tabl2.str[5,1]),]
    
    est[4,i] <- with(temp.rd, mean(outcome.1, na.rm=T))
    est[5,i] <- with(temp.rd, sd(outcome.1, na.rm=T))
    est[6,i] <- with(temp.rd, sd(temp.rd$outcome.1[treat==0], na.rm=T))
    
    est[7,i] <- nrow(temp.rd)
    est[8,i] <- with(temp.rd, table(non.zero))[2]
    
    two.sided.pvalue[i,1] <- robust$pv[3]
    one.sided.pvalue[i,1] <- pnorm(abs(robust$z[3]), lower.tail=FALSE)   
  }
   
  colnames(est) <- c("Local Poly Robust, 2000", "Local Poly Robust, 2004", 
                     "Local Poly Robust, 2008", "Local Poly Robust, 2012")
  rownames(est) <- c("Estimates", 
                     "Std. error", 
                     "t-stat",
                     "Mean dep. var.", 
                     "Standard-deviation dep. var.",
                     "Standard-deviation dep. var. (control)",
                     "n", "non zero") 
  table <- xtable(est,
                  align="r|cccc", 
                  include.rownames=TRUE)
  return(list(table, two.sided.pvalue, one.sided.pvalue)) 
}

############## Conjoint plot

plot_cjoint <- function(results=results, x_axis = c(-.5, .5)){
  temp <- summary(results)
  #adjust <- .15
  y.axis <- c(1, 2, 3, 4, 5, 6)
  coef.vec <- c(0, temp$amce[1,3], 0, temp$amce[2,3], temp$amce[3,3], temp$amce[4,3])
  
  lower.bounds <- c(0, temp$amce[1,3]-(qnorm(0.975)*temp$amce[1,4]),
                    0, temp$amce[2,3]-(qnorm(0.975)*temp$amce[2,4]), 
                    temp$amce[3,3]-(qnorm(0.975)*temp$amce[3,4]), 
                    temp$amce[4,3]-(qnorm(0.975)*temp$amce[4,4]))
  upper.bounds <- c(0, temp$amce[1,3]+(qnorm(0.975)*temp$amce[1,4]),
                    0, temp$amce[2,3]+(qnorm(0.975)*temp$amce[2,4]), 
                    temp$amce[3,3]+(qnorm(0.975)*temp$amce[3,4]), 
                    temp$amce[4,3]+(qnorm(0.975)*temp$amce[4,4]))
  par(mar=c(4, 9, 2, 0.5))
  
  plot(coef.vec, y.axis, type="p", axes = F, xlab = "", ylab = "", 
       pch = 16, cex = 1,
       xlim= x_axis, col=ifelse(y.axis < 3,'gray65','black'),
       main = "", cex.main=1) 
  abline(v=0, lty=3)
  #+adjust
  segments(lower.bounds, y.axis, upper.bounds, y.axis, lwd =  2, col=ifelse(y.axis < 3,'gray30','black')) #gray65, gray8
  points(coef.vec, 1:length(coef.vec), pch = 16, cex = 1, col=ifelse(y.axis < 3,'gray30','black')) 
  
  names <- c("Reference:\nNo Joint Provision", "Yes Joint Provision", 
             "Reference:\nState Provider", "State Provider,\nFed. Funds", "Non-State Provider", 
             "Non-State Provider,\nFed. Funds")
  axis(2, at = seq(1, 6), labels=names, 
       las = 1)
  axis(1, at = seq(x_axis[1], x_axis[2],  by=.25),  labels= seq(x_axis[1], x_axis[2],  by=.25))
}

#############log transformation

inv_hyp <- function(y){
  re <- log(y+(((y^2)+1)^(1/2)))
  return(re)}

###########  Rdd Table Balance
rdd.table.balance <- function(runvar=NULL, treat=NULL,
                              outcome= NULL, 
                              title="Covariate"){
  
  all_covariates <- list()
  for (j in 1:length(outcome)){
    
    dep_var <- outcome[[j]] 
    temp <- as.data.frame(cbind(runvar, treat, dep_var))
    
    print(title[j])
    bw <- c(0.005, 0.01, 0.02, 0.03, 0.04, 0.05)    
    cell <- matrix(NA, 4, length(bw))
    
    for (i in 1:length(bw)){
      temp1 <- temp[abs(temp$runvar)<=bw[i],]
      
      model <- with(temp1, lm(dep_var ~ treat))
      
      cell[1:2,i] <- c(summary(model)$coefficients[2,1], 
                       sqrt(diag(vcovHC(model, type = "HC1")))[2])
      
      Ys <- with(temp1, genouts(dep_var, treat, ate=0))
      probs <- with(temp1, genprobexact(treat))
      perms <- with(temp1, genperms(treat, maxiter=100000))
      ate <- with(temp1, estate(dep_var, treat, prob=probs))
      distout <- with(temp1, gendist(Ys, perms=perms, prob=probs))
      resultp <- dispdist(distout, ate)
      
      cell[3,i] <- resultp$two.tailed.p.value.abs   
      cell[4,i] <- ifelse(ate > 0, resultp$greater.p.value, resultp$lesser.p.value)
      rownames(cell) <- c(title[j], paste("standard error", title[j]), 
                          paste("two-tailed p-value", title[j]),
                          paste("one-tailed p-value", title[j]))  
    }
    all_covariates[[j]] <- cell
  }
  return(all_covariates) 
}  

### Local Polynomial Balance

rdd.table.balance.poly <- function(runvar=NULL, treat=NULL,
                                    outcome= NULL, 
                                    title="Covariate"){
  all_covariates <- list()
  for (j in 1:length(outcome)){
    
    dep_var <- outcome[[j]] 
    cell <- summary(rdrobust(x = runvar, y = dep_var, c = 0))$coefficients[3,]
    all_covariates[[j]] <- cell
  }
  return(all_covariates)
}


################################################################################
# RDDENSITY R PACKAGE -- rddensity -- manipulation test
# Authors: Matias D. Cattaneo, Michael Jansson, Xinwei Ma
################################################################################
#!version 0.1  31Mar2015

S.generate <- function(p, low=-1, up=1, kernel="triangular") {
  popwarning <- FALSE
  S <- matrix(rep(0, (p+1)^2), ncol=(p+1))
  for (i in 1:(p+1)) {
    for (j in 1:(p+1)) {
      if (kernel == "uniform") { 
        integrand <- function(x) { x^(i+j-2)*0.5 } 
      } else if (kernel == "epanechnikov") {
        integrand <- function(x) { x^(i+j-2)*0.75*(1-x^2) } 
      } else {
        integrand <- function(x) { x^(i+j-2)*(1-abs(x)) } 
      }  
      S[i,j] <- (integrate(integrand, lower=low, upper=up)$value)                             
    }
  }
  if (popwarning) {warning(text)}
  return(S)
}

S.minus.generate <- function(p, kernel="triangular") {
  S <- S.generate(p, low=-1, up=0, kernel=kernel)
  temp <- matrix(0, ncol=p+2, nrow=p+2)
  temp[1:2, 1:2] <- S[1:2, 1:2]
  if (p>1) {
    temp[1:2, 4:(p+2)] <- S[1:2, 3:(p+1)]
    temp[4:(p+2), 1:2] <- S[3:(p+1), 1:2]
    temp[4:(p+2), 4:(p+2)] <- S[3:(p+1), 3:(p+1)]
  }
  return(temp)
}

S.plus.generate <- function(p, kernel="triangular") {
  S <- S.generate(p, low=0, up=1, kernel=kernel)
  temp <- matrix(0, ncol=p+2, nrow=p+2)
  temp[1, 1] <- S[1, 1]
  temp[1, 3:(p+2)] <- S[1, 2:(p+1)]
  temp[3:(p+2), 1] <- S[2:(p+1), 1]
  temp[3:(p+2), 3:(p+2)] <- S[2:(p+1), 2:(p+1)]
  return(temp)
}

C.generate <- function(k, p, low=-1, up=1, kernel="triangular") {
  popwarning <- FALSE
  C <- matrix(rep(0, (p+1)), ncol=1)
  for (i in 1:(p+1)) {
    if (kernel == "uniform") { 
      integrand <- function(x) { x^(i+k-1)*0.5 } 
    } else if (kernel == "epanechnikov") {
      integrand <- function(x) { x^(i+k-1)*0.75*(1-x^2) }
    }
    else { 
      integrand <- function(x) { x^(i+k-1)*(1-abs(x)) } 
    }
    C[i,1] <- (integrate(integrand, lower=low, upper=up)$value)                             
  }
  if (popwarning) {warning(text)}
  return(C)
}

C.minus.generate <- function(k, p, kernel="triangular") {
  C <- C.generate(k=k, p=p, kernel=kernel, low=-1, up=0)
  temp <- matrix(0, ncol=1, nrow=p+2)
  temp[1:2] <- C[1:2]
  if(p>1) {
    temp[4:(p+2)] <- C[3:(p+1)]
  }
  return (temp)
}

C.plus.generate <- function(k, p, kernel="triangular") {
  C <- C.generate(k=k, p=p, kernel=kernel, low=0, up=1)
  temp <- matrix(0, ncol=1, nrow=p+2)
  temp[1] <- C[1]
  temp[3:(p+2)] <- C[2:(p+1)]
  return (temp)
}

G.generate <- function(p, low=-1, up=1, kernel="triangular") {
  popwarning <- FALSE
  G <- matrix(rep(0, (p+1)^2), ncol=(p+1))
  for (i in 1:(p+1)) {
    for (j in 1:(p+1)) {
      if (kernel == "uniform") {
        G[i,j] <- integrate(function(y) { 
          sapply(y, function(y) {
            integrate(function(x) x^i * y^(j-1)*0.25, low, y)$value
          })
        }, low, up)$value + 
          integrate(function(y) { 
            sapply(y, function(y) {
              integrate(function(x) x^(i-1) * y^j*0.25, y, up)$value
            })
          }, low, up)$value
      } else if (kernel == "epanechnikov") {
        G[i,j] <- integrate(function(y) { 
          sapply(y, function(y) {
            integrate(function(x) x^i * y^(j-1) * 0.75^2 *
                        (1-x^2) * (1-y^2), low, y)$value
          })
        }, low, up)$value + 
          integrate(function(y) { 
            sapply(y, function(y) {
              integrate(function(x) x^(i-1) * y^j * 0.75^2 *
                          (1-x^2) * (1-y^2), y, up)$value
            })
          }, low, up)$value
      } else {
        G[i,j] <- integrate(function(y) { 
          sapply(y, function(y) {
            integrate(function(x) x^i * y^(j-1) *
                        (1-abs(x)) * (1-abs(y)), low, y)$value
          })
        }, low, up)$value + 
          integrate(function(y) { 
            sapply(y, function(y) {
              integrate(function(x) x^(i-1) * y^j * 
                          (1-abs(x)) * (1-abs(y)), y, up)$value
            })
          }, low, up)$value
      }    
    }
  }
  if (popwarning) {warning(text)}
  return(G)
}

G.minus.generate <- function(p, kernel="triangular") {
  G <- G.generate(p, low=-1, up=0, kernel=kernel)
  temp <- matrix(0, ncol=p+2, nrow=p+2)
  temp[1:2, 1:2] <- G[1:2, 1:2]
  if (p>1) {
    temp[1:2, 4:(p+2)] <- G[1:2, 3:(p+1)]
    temp[4:(p+2), 1:2] <- G[3:(p+1), 1:2]
    temp[4:(p+2), 4:(p+2)] <- G[3:(p+1), 3:(p+1)]
  }
  return(temp)
}

G.plus.generate <- function(p, kernel="triangular") {
  G <- G.generate(p, low=0, up=1, kernel=kernel)
  temp <- matrix(0, ncol=p+2, nrow=p+2)
  temp[1, 1] <- G[1, 1]
  temp[1, 3:(p+2)] <- G[1, 2:(p+1)]
  temp[3:(p+2), 1] <- G[2:(p+1), 1]
  temp[3:(p+2), 3:(p+2)] <- G[2:(p+1), 2:(p+1)]
  return(temp)
}

Psi.generate <- function(p) {
  if (p > 1) {
    temp <- c(1, 0, 0, (-1)^(2:p))
  } else {
    temp <- c(1, 0, 0)
  }
  temp <- diag(temp)
  temp[2, 3] <- temp[3, 2] <- -1
  return(temp)
}

################################################################################
# RDDENSITY R PACKAGE -- rddensity -- Mata functions
# Authors: Matias D. Cattaneo, Michael Jansson, Xinwei Ma
################################################################################
#!version 0.1  31Mar2015

## NOTE: DATA IS ASSUMED TO BE IN ASCENDING ORDER

rddensity_fV <- function(Y, X, Nl, Nr, Nlh, Nrh, hl, hr, p, s, kernel, fitselect) {
  
  N <- Nl + Nr; Nh <- Nlh + Nrh
  Y <- matrix(Y, ncol=1); X <- matrix(X, ncol=1)
  W <- rep(NA, Nh)
  if (kernel == "uniform") {
    W[1:Nlh] <- 1 / (2 * hl); W[(Nlh+1):Nh] <- 1 / (2 * hr)
  } else if (kernel == "triangular") {
    W[1:Nlh] <- (1 + X[1:Nlh]/hl) / hl; W[(Nlh+1):Nh] <- (1 - X[(Nlh+1):Nh]/hr) / hr
  } else {
    W[1:Nlh] <- 0.75 * (1 - (X[1:Nlh]/hl)^2) / hl; W[(Nlh+1):Nh] <- 0.75 * (1 - (X[(Nlh+1):Nh]/hr)^2) / hr
  }
  
  if (fitselect == "restricted") {
    Xp <- matrix(NA, ncol=p+2, nrow=Nh)
    Xp[, 1] <- 1
    Xp[1:Nlh, 2] <- X[1:Nlh] / hl; Xp[(Nlh+1):Nh, 2] <- 0
    Xp[1:Nlh, 3] <- 0; Xp[(Nlh+1):Nh, 3] <- X[(Nlh+1):Nh] / hr
    if (p>1) {
      for (j in 4:(p+2)) {
        Xp[1:Nlh, j] <- (X[1:Nlh] / hl)^(j-2); Xp[(Nlh+1):Nh, j] <- (X[(Nlh+1):Nh] / hr)^(j-2)
      }
      v <- c(0, 1, 1, 2:p)
    } else {
      v <- c(0, 1, 1)
    }
    Hp <- diag(hl^v)
  } else {
    Xp <- matrix(NA, ncol=2*p+2, nrow=Nh)
    Hp <- rep(NA, 2*p+2)
    for (j in 1:(2*p+2)) {
      if (j %% 2) {
        Xp[1:Nlh, j] <- (X[1:Nlh] / hl)^((j-1)/2); Xp[(Nlh+1):Nh, j] <- 0
        Hp[j] <- hl^((j-1)/2)
      } else {
        Xp[1:Nlh, j] <- 0; Xp[(Nlh+1):Nh, j] <- (X[(Nlh+1):Nh] / hr)^((j-2)/2)
        Hp[j] <- hr^((j-2)/2)
      }
    }
    Hp <- diag(Hp)
  }
  
  out <- matrix(NA, ncol=4, nrow=4)
  colnames(out) <- c("hat", "jackknife", "plugin", "s")
  rownames(out) <- c("l", "r", "diff", "sum")
  
  Sh <- t(Xp) %*% diag(W) %*% Xp 
  Sinv <- try(solve(Sh), silent=TRUE)
  if (typeof(Sinv) == "character") {
    return(data.frame(out))
  }
  
  XpWY <- t(Xp) %*% diag(W) %*% Y
  b <- solve(Hp) %*% Sinv %*% XpWY
  
  if (fitselect == "restricted") {
    out[1, 1] <- b[2]; out[2, 1] <- b[3]; out[3, 1] <- b[3] - b[2]; out[4, 1] <- b[3] + b[2];
    out[1, 4] <- out[2, 4] <- b[s+2]; out[3, 4] <- 0; out[4, 4] <- 2 * out[1, 4] 
  } else {
    out[1, 1] <- b[3]; out[2, 1] <- b[4]; out[3, 1] <- b[4] - b[3]; out[4, 1] <- b[4] + b[3];
    out[1, 4] <- b[2*s+1]; out[2, 4] <- b[2*s+2]; out[3, 4] <- out[2, 4] - out[1, 4]; out[4, 4] <- out[2, 4] + out[1, 4]
  }
  
  # Jackknife
  XpW <- diag(W) %*% Xp; L <- matrix(NA, nrow=dim(Xp)[1], ncol=dim(Xp)[2])
  L[1, ] <- colSums(XpW[2:Nh, ]) / (N - 1)
  for (i in 2:Nh) {
    L[i, ] <- L[i-1, ] - XpW[i, ] / (N - 1)
  }
  V = solve(Hp) %*% Sinv %*% (t(L) %*% L) %*% Sinv %*% solve(Hp) 
  if (fitselect == "restricted") {
    out[1, 2] <- V[2, 2]; out[2, 2] <- V[3, 3]; out[3, 2] <- V[2, 2] + V[3, 3] - 2 * V[2,3]; out[4, 2] <- V[2, 2] + V[3, 3] + 2 * V[2,3]
  } else {
    out[1, 2] <- V[3, 3]; out[2, 2] <- V[4, 4]; out[3, 2] <- V[3, 3] + V[4, 4] - 2 * V[3,4]; out[4, 2] <- V[3, 3] + V[4, 4] + 2 * V[3,4]
  }
  
  # plugin
  if (fitselect=="unrestricted") {
    S <- S.generate(p, low=0, up=1, kernel=kernel)
    G <- G.generate(p, low=0, up=1, kernel=kernel)
    V <- solve(S) %*% G %*% solve(S)
    out[1, 3] <- out[1, 1] * V[2, 2] / (N * hl);  out[2, 3] <- out[2, 1] * V[2, 2] / (N * hr); out[3, 3] <- out[4, 3] <- out[1, 3] + out[2, 3]
  } else {
    S <- S.plus.generate(p=p, kernel=kernel)
    G <- G.plus.generate(p=p, kernel=kernel)
    Psi <- Psi.generate(p=p)
    Sm <- Psi %*% S %*% Psi; Gm <- Psi %*% G %*% Psi 
    V <- solve(out[1, 1] * Sm + out[2, 1] * S) %*% (out[1, 1]^3 * Gm + out[2, 1]^3 * G) %*% solve(out[1, 1] * Sm + out[2, 1] * S)
    out[1, 3] <- V[2, 2] / (N * hl); out[2, 3] <- V[3, 3] / (N * hl); out[3, 3] <- (V[2, 2] + V[3, 3] - 2 * V[2,3]) / (N * hl); out[4, 3] <- (V[2, 2] + V[3, 3] + 2 * V[2,3]) / (N * hl)
  }
  
  for (i in 1:4) {
    for (j in 2:3) {
      if (out[i, j] < 0) { out[i,j] <- NA }
    }
  }
  
  return(data.frame(out))
}

h.opt.density <- function(x, p, N, dgp_F1, dgp_Fp1, f_low, f_up, kernel="triangular") {
  
  # warning message
  if ((kernel != "triangular") & (kernel != "epanechnikov") & (kernel != "uniform")) {
    text <- paste("No kernel as ", toString(kernel), 
                  ", triangular kernel (default) is used.", sep="")
    message(text)
    kernel <- "triangular"
  }
  # end of warning message
  
  if (x==f_low | x==f_up) {
    c_low <- 0
  } else {
    c_low <- -1
  }
  c_up <- 1
  
  e <- matrix(rep(0, p+1), ncol=1)
  e[2] <- 1
  
  S1    <- S.generate(p=p, low=c_low, up=c_up, kernel=kernel)
  Cp1   <- C.generate(k=p+1, p=p, low=c_low, up=c_up, kernel=kernel)
  G     <- G.generate(p=p, low=c_low, up=c_up, kernel=kernel)
  
  kappa <- N^(-1/(2*p+1)) * (t(e) %*% solve(S1) %*% G %*% solve(S1) %*% e)^(1/(2*p+1)) * 
    (abs(t(e) %*% solve(S1) %*% Cp1))^(-2/(2*p+1)) * factorial(p+1)^(2/(2*p+1)) * (2*p)^(-1/(2*p+1))
  
  biassq <- (abs(dgp_Fp1))^(-2/(2*p+1))
  var1 <-  (dgp_F1)^(1/(2*p+1)) 
  
  h.opt <- kappa * biassq * var1
  return(as.numeric(h.opt))
}

h.opt.density.res <- function(p, N, dgp_F1_l, dgp_F1_r, dgp_Fp1_l, dgp_Fp1_r, kernel="triangular") {
  
  # warning message
  if ((kernel != "triangular") & (kernel != "epanechnikov") & (kernel != "uniform")) {
    text <- paste("No kernel as ", toString(kernel), 
                  ", triangular kernel (default) is used.", sep="")
    message(text)
    kernel <- "triangular"
  }
  # end of warning message
  
  e <- matrix(rep(0, p+2), ncol=1)
  if (p %% 2) {
    e[2] <- -1; e[3] <- 1
  } else {
    e[2] <- 1; e[3] <- 1
  }
  
  S.minus <- S.minus.generate(p=p, kernel=kernel); S.plus <- S.plus.generate(p=p, kernel=kernel)
  S.inv    <- solve(S.minus * dgp_F1_l + S.plus * dgp_F1_r)
  C   <- dgp_F1_l * dgp_Fp1_l * C.minus.generate(k=p+1, p=p, kernel=kernel) +
    dgp_F1_r * dgp_Fp1_r * C.plus.generate(k=p+1, p=p, kernel=kernel)
  G       <- dgp_F1_l^3 * G.minus.generate(p=p, kernel=kernel) + dgp_F1_r^3 * G.plus.generate(p=p, kernel=kernel) + 
    dgp_F1_l^2 * dgp_F1_r * tcrossprod(S.minus[2, ], S.plus[1, ]) +
    dgp_F1_l^2 * dgp_F1_r * tcrossprod(S.plus[1, ], S.minus[2, ])
  
  h.opt <- N^(-1/(2*p+1)) * (t(e) %*% S.inv %*% G %*% S.inv %*% e)^(1/(2*p+1)) * 
    (abs(t(e) %*% S.inv %*% C))^(-2/(2*p+1)) * factorial(p+1)^(2/(2*p+1)) * (2*p)^(-1/(2*p+1))
  
  return(as.numeric(h.opt))
}

rddensity_H <- function(x, p){
  if (p==0)  out = 1
  if (p==1)  out = x
  if (p==2)  out = x^2 - 1
  if (p==3)  out = x^3 - 3*x
  if (p==4)  out = x^4 - 6*x^2 + 3
  if (p==5)  out = x^5 - 10*x^3 + 15*x
  if (p==6)  out = x^6 - 15*x^4 + 45*x^2 - 15
  if (p==7)  out = x^7 - 21*x^5 + 105*x^3 - 105*x
  if (p==8)  out = x^8 - 28*x^6 + 210*x^4 - 420*x^2 + 105
  if (p==9)  out = x^9 - 36*x^7 + 378*x^5 - 1260*x^3 + 945*x
  if (p==10) out = x^10 - 45*x^8 + 630*x^6 - 3150*x^4 + 4725*x^2 - 945
  return(out)
}

#!version 0.1  31Mar2015

#source(paste0(dir, "replication/codes/rdbwdensity.R"))

rddensity <- function(X, c=0, p=2, q=0, kernel="", fitselect="", hl=0, hr=0, hscale=1, bwselect="", vce="", print.screen=TRUE) {
  
  ################################################################################
  # default values 
  ################################################################################
  if (q==0) { q <- p+1 }
  if (kernel == "") { kernel <- "triangular" }
  kernel <- tolower(kernel)
  if (fitselect == "") { fitselect <- "unrestricted" }
  fitselect <- tolower(fitselect)
  if (bwselect == "") { bwselect <- "comb" }
  bwselect <- tolower(bwselect)
  if (vce == "") { vce <- "jackknife" }
  vce <- tolower(vce)
  # end of default values
  
  ################################################################################
  # sample sizes
  ################################################################################
  X <- sort(X)
  N <- length(X); Nl <- sum(X<c); Nr <- sum(X>=c); Xmin <- min(X); Xmax <- max(X)
  # end of sample sizes
  
  ################################################################################
  # error handling
  ################################################################################
  if (c <= Xmin | c >= Xmax) { stop("The cutoff should be set within the range of the data.") }
  if (Nl <= 10 | Nr <= 10) { stop("Not enough observations to perform calculations.") }
  if (p!=1 & p!=2 & p!=3 & p!=4 & p!=5 & p!=6 & p!= 7) { stop("p must be an integer between 1 and 7.") }
  if (p > q) { stop("p cannot exceed q.") }
  if (kernel!="uniform" & kernel!= "triangular" & kernel!="epanechnikov") { stop("kernel incorrectly specified.") }
  if (fitselect!="unrestricted" & fitselect!="restricted") { stop("fitselect incorrectly specified.") }
  if (hl<0 | hr<0) { stop("Bandwidth cannot be negative.") }
  if (fitselect=="restricted" & hl!=hr) { stop("Bandwidths must be equal in restricted model.") }
  if (hscale<=0 | hscale>1) { stop("hscale must be positive and no larger than 1.") }
  if (bwselect!="each" & bwselect!="diff" & bwselect!="sum" & bwselect!="comb") { stop("bwselect incorrectly specified.") }
  if (fitselect=="restricted" & bwselect=="each") { stop("bwselect=each is not available in the restricted model.") }
  if (vce!="plugin" & vce!="jackknife") { stop("vce incorrectly specified.") }
  # end of error handling
  
  ################################################################################
  # bandwidth selection
  ################################################################################
  if (hl > 0) { bwselectl <- "mannual" } else { bwselectl <- "estimated" }
  if (hr > 0) { bwselectr <- "mannual" } else { bwselectr <- "estimated" }
  if (hl==0 | hr==0) {
    out <- rdbwdensity(X=X, c=c, p=p, kernel=kernel, fitselect=fitselect, vce=vce, print.screen=FALSE)
    if (fitselect=="unrestricted" & bwselect=="each" & hl==0)  hl = out[1,1]
  	if (fitselect=="unrestricted" & bwselect=="each" & hr==0)  hr = out[2,1]
		if (fitselect=="unrestricted" & bwselect=="diff" & hl==0)  hl = out[3,1]
		if (fitselect=="unrestricted" & bwselect=="diff" & hr==0)  hr = out[3,1]
		if (fitselect=="unrestricted" & bwselect=="sum"  & hl==0)  hl = out[4,1]
		if (fitselect=="unrestricted" & bwselect=="sum"  & hr==0)  hr = out[4,1]
		if (fitselect=="unrestricted" & bwselect=="comb" & hl==0)  hl = median(c(out[1,1], out[3,1], out[4,1]))
		if (fitselect=="unrestricted" & bwselect=="comb" & hr==0)  hr = median(c(out[2,1], out[3,1], out[4,1]))

		if (fitselect=="restricted" & bwselect=="diff" & hl==0)  hl = out[3,1]
		if (fitselect=="restricted" & bwselect=="diff" & hr==0)  hr = out[3,1]
		if (fitselect=="restricted" & bwselect=="sum"  & hl==0)  hl = out[4,1]
		if (fitselect=="restricted" & bwselect=="sum"  & hr==0)  hr = out[4,1]
		if (fitselect=="restricted" & bwselect=="comb" & hl==0)  hl = min(c(out[3,1], out[4,1]))
		if (fitselect=="restricted" & bwselect=="comb" & hr==0)  hr = min(c(out[3,1], out[4,1]))
  }
  # end of bandwidth selection
  
  ################################################################################
  # data trimming
  ################################################################################
  X <- X - c; Y <- (0:(N-1)) / (N-1); 
  Xh <- X[(X>=-1*hl*hscale) & (X<=hr*hscale)]; Yh <- Y[(X>=-1*hl*hscale) & (X<=hr*hscale)]
  Nlh <- sum(Xh < 0); Nrh <- sum(Xh >= 0)
  if (Nlh < 5 | Nrh < 5) { stop("Not enough observations to perform calculation.") }
  Nh <- Nlh + Nrh
  # end of data trimming
  
  ################################################################################
  # estimation
  ################################################################################
  fV_q <- rddensity_fV(Y=Yh, X=Xh, Nl=Nl, Nr=Nr, Nlh=Nlh, Nrh=Nrh, hl=hl*hscale, hr=hr*hscale, p=q, s=1, kernel=kernel, fitselect=fitselect)
  T_asy <- fV_q[3, 1] / sqrt(fV_q[3, 3]); T_jk <- fV_q[3, 1] / sqrt(fV_q[3, 2])
  p_asy <- pnorm(abs(T_asy), lower.tail=FALSE) * 2; p_jk <- pnorm(abs(T_jk), lower.tail=FALSE) * 2
  
  result <- list( hat=   list(left=fV_q[1,1], right=fV_q[2,1], diff=fV_q[3,1]), 
                  sd.asy=list(left=sqrt(fV_q[1,3]), right=sqrt(fV_q[2,3]), diff=sqrt(fV_q[3,3])), 
                  sd.jk= list(left=sqrt(fV_q[1,2]), right=sqrt(fV_q[2,2]), diff=sqrt(fV_q[3,2])),
                  test=  list(t.asy=T_asy, t.jk=T_jk, p.asy=p_asy, p.jk=p_jk), 
                  N=     list(full=N, left=Nl, right=Nr, eff.left=Nlh, eff.right=Nrh), 
                  h=     list(left=hl*hscale, right=hr*hscale))
  
  if (print.screen) {
    cat("\nRD Manipulation Test using local polynomial density estimation.\n")
    cat("\n")
    
    cat(paste(format("Number of obs =", width=30), toString(N), sep="")); cat("\n")
    cat(paste(format("Model =", width=30), fitselect, sep="")); cat("\n")
    cat(paste(format("Kernel =", width=30), kernel, sep="")); cat("\n")
    if (bwselectl!="mannual" | bwselectr!="mannual") {
      cat(paste(format("BW method =", width=30), bwselect, sep="")); cat("\n")
    }
    if (hscale < 1) {
      cat(paste(format("BW scale =", width=30), toString(round(hscale, 3)), sep="")); cat("\n")
    }
    cat(paste(format("VCE method =", width=30), vce, sep="")); cat("\n")
    cat("\n")
    
    cat(paste(format(paste("Cutoff c = ", toString(round(c, 3)), sep=""), width=30), format("Left of c", width=20), format("Right of c", width=20), sep="")); cat("\n")
    cat(paste(format("Number of obs", width=30), format(toString(Nl), width=20), format(toString(Nr), width=20), sep="")); cat("\n")
    cat(paste(format("Eff. Number of obs", width=30), format(toString(Nlh), width=20), format(toString(Nrh), width=20), sep="")); cat("\n")
    cat(paste(format("Min Running var.", width=30), format(toString(round(min(X[X<0]+c), 3)), width=20), format(toString(round(min(X[X>=0]+c), 3)), width=20), sep="")); cat("\n")
    cat(paste(format("Max Running var.", width=30), format(toString(round(max(X[X<0]+c), 3)), width=20), format(toString(round(max(X[X>=0]+c), 3)), width=20), sep="")); cat("\n")
    cat(paste(format("Order loc. poly. (p)", width=30), format(toString(p), width=20), format(toString(p), width=20), sep="")); cat("\n")
    if (q > p) {
      cat(paste(format("Order BC (q)", width=30), format(toString(q), width=20), format(toString(q), width=20), sep="")); cat("\n")
    } else {
      cat(paste(format("No bias correction (q)", width=30), format(toString(q), width=20), format(toString(q), width=20), sep="")); cat("\n")
    }
    
    cat(paste(format("Bandwidths (hl,hr)", width=30), format(bwselectl, width=20), format(bwselectr, width=20), sep="")); cat("\n")
    if (hscale < 1) {
      cat(paste(format("Bandwidth values", width=30), format(paste(toString(round(hl, 3)), " * ", toString(round(hscale, 3)), sep=""), width=20), format(paste(toString(round(hr, 3)), " * ", toString(round(hscale, 3)), sep=""), width=20), sep="")); cat("\n")
    } else {
      cat(paste(format("Bandwidth values", width=30), format(toString(round(hl, 3)), width=20), format(toString(round(hr, 3)), width=20), sep="")); cat("\n")
    }
    cat("\n")
    
    cat(paste(format("Method", width=30), format("T", width=20), format("P > |T|", width=20), sep="")); cat("\n")
    if (p == q & hscale==1) {
      if (vce == "plugin") {
        cat(paste(format("Conventional", width=30), format(toString(round(T_asy, 4)), width=20), format(toString(round(p_asy, 4)), width=20), sep="")); cat("\n")
      } else {
        cat(paste(format("Conventional", width=30), format(toString(round(T_jk, 4)), width=20), format(toString(round(p_jk, 4)), width=20), sep="")); cat("\n")
      }      
    } else if (p == q & hscale<1) {
      if (vce == "plugin") {
        cat(paste(format("Undersmoothed", width=30), format(toString(round(T_asy, 4)), width=20), format(toString(round(p_asy, 4)), width=20), sep="")); cat("\n")
      } else {
        cat(paste(format("Undersmoothed", width=30), format(toString(round(T_jk, 4)), width=20), format(toString(round(p_jk, 4)), width=20), sep="")); cat("\n")
      }  
    } else if (p < q & hscale == 1) {
      if (vce == "plugin") {
        cat(paste(format("Robust Bias-Corrected", width=30), format(toString(round(T_asy, 4)), width=20), format(toString(round(p_asy, 4)), width=20), sep="")); cat("\n")
      } else {
        cat(paste(format("Robust Bias-Corrected", width=30), format(toString(round(T_jk, 4)), width=20), format(toString(round(p_jk, 4)), width=20), sep="")); cat("\n")
      }  
    } else {
      if (vce == "plugin") {
        cat(paste(format("Robust BC, Undersmoothed", width=30), format(toString(round(T_asy, 4)), width=20), format(toString(round(p_asy, 4)), width=20), sep="")); cat("\n")
      } else {
        cat(paste(format("Robust BC, Undersmoothed", width=30), format(toString(round(T_jk, 4)), width=20), format(toString(round(p_jk, 4)), width=20), sep="")); cat("\n")
      }  
    }
    cat("\n")
    if (p==q & hscale==1) {
      warning("Without bias correction or undersmoothing, the size of the test may not be correct.") 
    }
    if (p==1 & q==1) {
      warning("With p=q=1, the size of the test may not be correct.")
    }
    cat("\n")
    
  } else {
    return(result)
  }
}



rdbwdensity <- function(X, c=0, p=2, kernel="", fitselect="", vce="", print.screen=TRUE) {
  
  ################################################################################
  # default values 
  ################################################################################
  if (kernel == "") { kernel <- "triangular" }
  kernel <- tolower(kernel)
  if (fitselect == "") { fitselect <- "unrestricted" }
  fitselect <- tolower(fitselect)
  if (vce == "") { vce <- "jackknife" }
  vce <- tolower(vce)
  # end of default values
  
  ################################################################################
  # sample sizes
  ################################################################################
  X <- sort(X)
  N <- length(X); Nl <- sum(X<c); Nr <- sum(X>=c); Xmin <- min(X); Xmax <- max(X)
  # end of sample sizes
  
  ################################################################################
  # error handling
  ################################################################################
  if (c <= Xmin | c >= Xmax) { stop("The cutoff should be set within the range of the data.") }
  if (Nl <= 10 | Nr <= 10) { stop("Not enough observations to perform calculations.") }
  if (p!=1 & p!=2 & p!=3 & p!=4 & p!=5 & p!=6 & p!= 7) { stop("p must be an integer between 1 and 7.") }
  if (kernel!="uniform" & kernel!= "triangular" & kernel!="epanechnikov") { stop("kernel incorrectly specified.") }
  if (fitselect!="unrestricted" & fitselect!="restricted") { stop("fitselect incorrectly specified.") }
  if (vce!="plugin" & vce!="jackknife") { stop("vce incorrectly specified.") }
  # end of error handling
  
  ################################################################################
  # select preliminary bandwidth
  ################################################################################
  X <- X - c; Xmu <- mean(X); Xsd <- sd(X)
  fhatb <- 1 / (rddensity_H(Xmu / Xsd, p+2)^2 * dnorm(Xmu / Xsd))
  fhatc <- 1 / (rddensity_H(Xmu / Xsd, p)^2 * dnorm(Xmu / Xsd))
  # these constants are for uniform kernel
  Cb <- c(25884.444444494150957,3430865.4551236177795,845007948.04262602329,330631733667.03808594,187774809656037.3125,145729502641999264,146013502974449876992)
  Cc <- c(4.8000000000000246914,548.57142857155463389,100800.00000020420703,29558225.458100609481,12896196859.612621307,7890871468221.609375,6467911284037581)
  bn <- ((2*p+1)/4 * fhatb * Cb[p] / N)^(1/(2*p+5))
  cn <- (1/(2*p) * fhatc * Cc[p] / N)^(1/(2*p+1))
  bn <- bn * Xsd; cn <- cn * Xsd
  # end of select preliminary bandwidth
  
  ################################################################################
  # estimate main bandwidth
  ################################################################################
  Y <- (0:(N-1)) / (N-1)
  Yb <- Y[abs(X) <= bn]; Xb <- X[abs(X) <= bn]; Yc <- Y[abs(X) <= cn]; Xc <- X[abs(X) <= cn]
  Nlb <- sum(Xb < 0); Nrb <- sum(Xb >= 0); Nlc <- sum(Xc < 0); Nrc <- sum(Xc >= 0)
  
  hn <- matrix(NA, ncol=3, nrow=4)
  colnames(hn) <- c("bw", "variance", "biassq"); rownames(hn) <- c("l", "r", "diff", "sum")
  fV_b <- rddensity_fV(Y=Yb, X=Xb, Nl=Nl, Nr=Nr, Nlh=Nlb, Nrh=Nrb, hl=bn, hr=bn, p=p+2, s=p+1, kernel=kernel, fitselect=fitselect)
  fV_c <- rddensity_fV(Y=Yc, X=Xc, Nl=Nl, Nr=Nr, Nlh=Nlc, Nrh=Nrc, hl=cn, hr=cn, p=p, s=1, kernel=kernel, fitselect=fitselect)
  
  if (vce == "plugin") { hn[, 2] <- N * cn * fV_c[, 3]  } else { hn[, 2] <- N * cn * fV_c[, 2] }
  if (fitselect == "unrestricted") {
    S <- S.generate(p=p, low=0, up=1, kernel=kernel); C <- C.generate(k=p+1, p=p, low=0, up=1, kernel=kernel)
    hn[1, 3] <- fV_b[1, 4] * (solve(S) %*% C)[2] * (-1)^p
    hn[2, 3] <- fV_b[2, 4] * (solve(S) %*% C)[2]
    hn[3, 3] <- hn[2, 3] - hn[1, 3]; hn[4, 3] <- hn[2, 3] + hn[1, 3]
  } else {
    Splus <- S.plus.generate(p=p, kernel=kernel); Cplus <- C.plus.generate(k=p+1, p=p, kernel=kernel); Psi <- Psi.generate(p=p)
    Sinv <- solve(fV_c[2, 1] * Splus + fV_c[1, 1] * Psi%*%Splus%*%Psi)
    C <- fV_b[1, 4] * (fV_c[2, 1] * Cplus + (-1)^(p+1) * fV_c[1, 1] * Psi%*%Cplus)
    temp <- Sinv%*%C
    hn[1, 3] <- temp[2]; hn[2, 3] <- temp[3]; hn[3, 3] <- hn[2, 3] - hn[1, 3]; hn[4, 3] <- hn[2, 3] + hn[1, 3]
  }
  
  hn[, 3] <- hn[, 3]^2
  hn[, 1] <- (1/(2*p) * hn[, 2] / hn[, 3] / N)^(1/(2*p+1))
  # end of estimate main bandwidth
  
  for (i in 1:4) {
    if (hn[i, 2] < 0) { hn[i, 1] <- 0; hn[i, 2] <- NA }
    if (is.na(hn[i, 1])) { hn[i, 1] <- 0 }
  }
  
  ################################################################################
  # output the results
  ################################################################################
  if (print.screen) {
    cat("\nBandwidth selection for manipulation testing.\n")
    cat("\n")
    
    cat(paste(format("Number of obs =", width=30), toString(N), sep="")); cat("\n")
    cat(paste(format("Model =", width=30), fitselect, sep="")); cat("\n")
    cat(paste(format("Kernel =", width=30), kernel, sep="")); cat("\n")
    cat(paste(format("VCE method =", width=30), vce, sep="")); cat("\n")
    cat("\n")
    
    cat(paste(format(paste("Cutoff c = ", toString(round(c, 3)), sep=""), width=30), format("Left of c", width=20), format("Right of c", width=20), sep="")); cat("\n")
    cat(paste(format("Number of obs", width=30), format(toString(Nl), width=20), format(toString(Nr), width=20), sep="")); cat("\n")
    cat(paste(format("Min Running var.", width=30), format(toString(round(min(X[X<0]+c), 3)), width=20), format(toString(round(min(X[X>=0]+c), 3)), width=20), sep="")); cat("\n")
    cat(paste(format("Max Running var.", width=30), format(toString(round(max(X[X<0]+c), 3)), width=20), format(toString(round(max(X[X>=0]+c), 3)), width=20), sep="")); cat("\n")
    cat(paste(format("Order loc. poly. (p)", width=30), format(toString(p), width=20), format(toString(p), width=20), sep="")); cat("\n")
    cat("\n")
    
    cat(paste(format("Target", width=30), format("Bandwidth", width=20), format("Variance", width=20), format("Bias^2", width=20), sep="")); cat("\n")
    cat(paste(format("left density", width=30), format(toString(round(hn[1,1], 4)), width=20), format(toString(round(hn[1,2], 4)), width=20), format(toString(round(hn[1,3], 4)), width=20), sep="")); cat("\n")
    cat(paste(format("right density", width=30), format(toString(round(hn[2,1], 4)), width=20), format(toString(round(hn[2,2], 4)), width=20), format(toString(round(hn[2,3], 4)), width=20), sep="")); cat("\n")
    cat(paste(format("difference densities", width=30), format(toString(round(hn[3,1], 4)), width=20), format(toString(round(hn[3,2], 4)), width=20), format(toString(round(hn[3,3], 4)), width=20), sep="")); cat("\n")
    cat(paste(format("sum densities", width=30), format(toString(round(hn[4,1], 4)), width=20), format(toString(round(hn[4,2], 4)), width=20), format(toString(round(hn[4,3], 4)), width=20), sep="")); cat("\n")
    cat("\n")
  } else {
    return(data.frame(hn))
  }
}














